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A  Note  On  The  Gaussian  Integration  Formula 

\ 

\  ABSTRACT 

The  Gaussian  integration  method  is  described  together  with  the  2-point 
implementation  for  numerical  integration.  A  FORTRAN  program  is  given  and  the 
result  of  a  test  example. 


ii 


BACKGROUND 

'i' 

1.  This  work  was  originally  completed  some  years  ago  and  the  software,  which 

is  incorporated  in  a  Maths  Library  on  the  UDT  computer  facility,  has  been  used 
extensively.  Gauss'  formula  provides  a  method  for  ^calculating  the  definite 
integral  of  a  given  function.  Au.  ’  r  . 

THEORY  '  "  ( 

2.  If  pm(x)i  pn(x)>  x  e  R,  are  polynomials  of  degree  m,n  respectively, 
with  m  >  n,  then  polynomials  q(x),  r(x)  exist  such  that: 

Pm(x)  -  Pn(x)q(x)  +  r(x)  1. 

where  the  degree  of  q(x)  is  m-n  and  the  degree  of  r(x)  is  n-1  at  most.  The 
polynomials  q(x),  r(x)  are  unique  and  called  the  quotient  and  remainder  of 
Pm(x)  with  reference  to  Pn(x). 

3.  Consider: 

+1 

I  P2n+l<x)dx 
-1 

and  let  pn+1(x)  be  a  Legendre  polynomial  of  degree  n+1.  Then  by  equation  1: 

+1  +1 

I  P2n+l(x)dx  "  I  {Pn+l(x)<1n(x)  +  rn(x)}dx  2‘ 

-1  -1 

But  any  polynomial  of  degree  n  can  be  expressed  as  a  linear  combination  of 
Legendre  polynomials  of  degree  k  <  n.  Thus 

n 

qn<*)  -  l  a|A(*)  3. 

k-0 

and  therefore,  using  the  orthogonal  properties  of  Legendre  functions: 

+1 

J  Pn+1(x)qn(x)dx  =0  4. 

-1 


gives: 


+  1 

1  P2n+1 


+1 

(x)dx  -  J  rn(x)dx 
-1 


5. 


1 


4. 


We  seek  solutions  of  the  form: 


+1  n 

I  P2n+l<x)dx  “  £  akP2n+l<xk> 
—1  k*0 


■  1  ak{Pn+l(xk)qn<xk)  +  rn(xk>} 
k=0 


but  equation  5  shows  the  integral  is  independant  of  the  choice  of  qn(x)  and 
therefore  equation  6  is  true  iff  xk  are  the  n+1  roots  of  Pn+^(x).  Let: 

n 

V(x)  =  TT  (x-xk>  7- 

k=0 


then  expanding  in  partial  fractions  gives: 


rn(x>  =  I  rn(xk> 

Vi(x)  k=0  “A+i(xk)(x"xk) 

Integration  of  equation  8  gives: 


8. 


J+rn(x)dx  -  J+1  E  Vl(x)rn(xk)dx  9. 

-1  -1  k'°  “n+l(xk)(x-xk> 

Therefore: 

+1  n 

J  P2n+1<x)dx  =  E  V2n*l<xk)  10- 

-1  k=0 

where: 

ak  =  f  Vl(x)dx  n. 

-1  Mn+l<xk)<x-xk) 

5.  If  we  use  the  Hermite  interpolation  formula  to  represent  a  function  f(x)  by 
a  polynomial  of  degree  2n+l,  the  error  in  the  integral  formula  is  given  by: 

en  '  f(2rU2)(x)  J  H  (x-xk)2dx  12. 

2n+2) !  -1  k-0 

where  -1  <  k  <  +1. 
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6.  There  are  n+1  values  of  x^,  and  hence  this  method  of  integration  is 
called  the  n+1  point  method.  For  the  Gaussian  2-point  integral  formula  put  n=l 
and  then  x^  k«0,l  are  the  zeros  of  P2(x),  ie  x^  are  given  by: 

(3x2-1)/2  -  0 


therefore: 


xk  =  ±1/43 


and: 


«n+i(x)  =  <x2  "  1/3) 
Hence  the  veighting  factors  are: 

+1 

a.  »  J  x  ±  1/43  dx 
-1  2x 

=  1 

Finally  the  formula  is: 


7.  In  general  f(x)dx  can  be  transformed  to  the  form  of  equation  13. 
Furthermore  the  accuracy  can  be  improved  by  subdividing  the  interval  a,b.  Let 
each  subinterval  have  half  vidth  h,  where: 

h  =  (b-a)/2n 


Then: 


b  n-1  2(k+l)h+a 

|  f(x)dx  «  £  J  f(x)dx 

a  k=.0  2kh+a 

n-1  +1 

=  Z  hj  f(x)dx' 
k-0  -1 

where: 

x  =■  x'h  +  h(2k+l)  +  a 

and  so  by  equation  13: 


3 


J  f(x)dx  »  h  l  {f(  a+h(2k+l-fl/')3))  +f(  a+h(2k+l-l/43))}+h.n.h4f(4)(ft)  14. 

a  k«0  135 

8.  The  error  term  In  equation  14  can  be  written: 

eR  -  (b-a)5f(4)(S) 

25n4135 


Taking  2n  intervals  the  error  will  be: 

e2n  =  (b-a)5£(4)(ii) 
2524n4135 
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Put: 


b 

J  £(x)dx  -  In  +  en 

a 

=  I2n  *  e2n 

then  a  closer  approximation  is  given  by: 
b 

I  f<x>dx  -  J2n  +  (I2n  '  V/15  15’ 

a 

9.  Gauss'  integral  formula  is  not  suitable  to  evaluate  an  integral  of  a 
function  given  in  tabular  form,  but  when  the  function  is  given  as  an  analytic 
expression  the  method  is  very  useful.  A  FORTRAN  procedure  is  described  to 
compute: 

b 

I  =  J  f(x)dx  +  e 
a 

where: 

I e|  S  tol 
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4-JUL-89  13:11 


qauss-int . for 


Pig-  1-01 


REAL  FUNCTION  GAUSS__INTEGRATION  (  F  r  A,  B  ,TOL  ) 

C  Dascrlption : 

c  Intagratas  th*  function  F(x)  ov»r  th*  ring*  (A,B),  whin  F  is  i 

C  rul  function  (to  bi  provided).  Thi  riputid  two-point  Gauss 

c  mathod  is  u» id  until  succassiv*  approxima t ions  diffir  by  liss 

c  thin  1 5 *TOL . 

c  Author: 

c  C  Richardson,  ARE  (Portland) 


C 

c 

c 


History : 

Issua  1.0  22  July  1902 

Mod  1.1  6  August  1985 


C  Modifications : 

c  1 . 1  To  adopt  coding  standard 


C 


c 


c 


Subroutln*  Arguaints: 

REAL  F 

A ,  - 
To  1 

Local  Variiblis: 

REAL  Pi , P2 , 

H, 

Cd ,  G 

INTEGER  Nc 

Main  Entry  Point: 

Nc-8 

DO  WHILE  ((Nc  .LT.  3 1 ) . OR . ( ABS ( P2-P1 ) 
P1-P2 

H« ( B— A ) / ( 2  *  Nc ) 

Cd-A-H 

G-H/l . 73205081 
P2-0. 

DO  M-l,Nc 
Cd-Cd+H+H 
P2-P2+F(Cd-G) 

P2*P2+F (  Cd+G ) 

ENDDO 

P2«P2*H 

Nc-2#Nc 

ENDDO 

GAUSS  INTEGRATION-P2+(P2-Pl )/15 . 

RETURN 

END 


l ( R )  Function 
MR)  Limits 

MR)  Raquirad  accuracy 


llntagration  astimataa 
tSubintarval  half  width 


iNumbar  of  subintarvals 


iNurabar  of  subintarvals 
.GT.  1 5 . *tol ) ) 

JSav*  last  astimata 
tSubintarval  half-width 


l Partial  sum 


lEstimatad  rasult 
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4-JUL-8  9  13:14 


GAUS 3 - I NT-TEST . F OR 


PROGRAM  GAUSS_INTEGRATIONJTEST 
Dtacr ipt ion : 

Test  Gaussian  integration  by  computing  the  integral  of  Sin(x) 
ovar  the  rang*  ( A , B ) . 


C  Richardson,  ARE  (Portland) 


Local  V*  r i a b 1 e e : 


EXTERNAL  T 


Main  Entry  Point: 


(Integration  Interval 
i Workspace 


TYPE  *,'The  integration  interval  A , B  must  be  specified' 
A* AS KIM  'Enter  value  of  A  (Degrees)') 

B«ASKR( 'Enter  value  of  B  (Degrees)') 


CALL  PLOT  INIT  ( ' WIND0W-A4-V ' )  Unitialise 

CALL  CHAR“SI2E(1.5,2.0) 

CALL  ORID“lIN_LIN( 30 . , 50.,180.,200.,A,-1.,B,1.,1) 

CALL  CURVE (  F  ,  A ,  8 , 1 . )  l Draw  funct 

CALL  MOVE ( ( A+B ) /2 . , — 1 . )  (Annotate  a 

CALL  CHAR  POSN ( -7 . 5 , - 2 . ) 

CALL  PLOT  TEXT( 'Angle  ( Deg r ee s ) ' , 1 5 ) 

CALL  MOVeTa,0.) 

CALL  CHAR_P0SN(-4 .  ,0  -  ) 

CALL  CHAR  ANGLE (90,) 

CALL  CHAR~POSN(-3 . ,0 . ) 

CALL  PLOT'"tEXT(  'Sin(X)  '  ,  6  ) 

CALL  CHAR^ANGLE ( 0 . ) 

P«GAUSS_INTEGRATION<F,A,8, .00001)  (Calculate 

P-P/57. 29578 

CALL  MOVE! A, l . )  (Annotate  g 

CALL  CHAR_POSN ( 0 . ,1.5) 

CALL  PLOT  TEXT ( 'Numeric  Integral- ', 1 8 ) 

CALL  PLOT~NUM(P, 3 ,4 ) 

P— COSD ( A ) — COSD ( B )  (Calculate 

CALL  MOVE ( A , 1 . )  (Annotate  g 

CALL  CHAR_POSN(0 . , 0 . 5 ) 

CALL  PLOT_TEXT( 'Analytic  Integral- 1 8 ) 

CALL  PLOT_NUM ( P , 3 , 4 ) 

CALL  OR1GIN(0 . , 0 . )  (Plot  title 

CALL  SCALE* 1.  ,1.  ) 

CALL  MOVE ( 30 . , 20 . ) 

CALL  CHARTS I ZE ( 1 . 8 , 2 . 4 ) 

CALL  PLOT_TRXT ( 

'Pig  1.  Gaussian  Integration  Of  The  Function  Sln(X)',51) 
ACCEPT  100, P 
CALL  PLOT  riN 
END 


(Initialise  graphics 


(Draw  function  over  A , B 
(Annotate  axes 


(Calculate  numerical  integral 
(Annotate  graph 


(Calculate  analytic  integral 
(Annotate  graph 


FUNCTION  F ( X ) 

F-SIND(X) 

RETURN 

END 


(A  typical  function 


